# DISPLAY MAPS WITH NEARMAP
# Note: most of this code is commented as it is only runnable with a NearMap API.

pacman::p_load(tidyverse, data.table, sf, cowplot, png, curl, viridis, qs, kableExtra)


## Preliminaries
rm(list = ls())

source("code/globals.R")

# NOT
# apikey<-readLines('D:\\jboomhower\\Dropbox\\Administrative\\Software\\Software\\NearMap\\NearMapAPIkey.txt',warn=FALSE)
#
# thisproj<-3310
#
#
#
# #------- Main homes dataset
# chars<-qread(file.path(WORKING,"maindata.qs")) %>%
#   select(ImportParcelID,BuildingOrImprovementNumber,incidentid,OUTSIDE,FIPS,
#          LotSizeAcres,destroyed,combinedyear,street,street100,street25,street100_byside,
#          incidentname)
#
# #-------- Spatial information
# homes<-qread(file.path(WORKING,"homepoints.qs")) %>%
#   rename(incidentid=UNIT_AND_NUMBER) %>%
#   merge(chars,by=c('ImportParcelID','BuildingOrImprovementNumber','incidentid')) %>%
#   st_transform(crs=4326) %>%
#   mutate(lon=st_coordinates(.)[,1],
#          lat=st_coordinates(.)[,2]) %>%
#   st_transform(thisproj)
# rm(chars)
#
# #--------- Buildings
# roofshapes<-qread(file.path(WORKING,'rooftop_locations.qs')) %>%
#   st_set_geometry('roofpolygon') %>%
#   dplyr::select(ImportParcelID,roofpolygon,roof_fail) %>%
#   st_transform(crs=st_crs(homes))
#
# #--------- Parcel boundary lines
# lotlines<-qread(file.path(WORKING,'rooftop_locations.qs')) %>%
#   st_set_geometry('geometry') %>%
#   dplyr::select(ImportParcelID,geometry) %>%
#   st_transform(crs=st_crs(homes))
#
# # FIre perimeters
# perims<-qread(file.path(WORKING,'cleanedperimeters.qs')) %>%
#   rename(incidentid=UNIT_AND_NUMBER) %>%
#   st_transform(crs=st_crs(homes))
#
#
# #------------------------------------------------------------------#
# #--- Helper functions to carry out steps of the figure making -----#
# #------------------------------------------------------------------#
#
# ############## GET MAP TILE NUMBERS #################################################
# #### Calculate Google Map Tile x and y values for a given lat,long point and zoom
# #https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Lon..2Flat._to_tile_numbers_2
#
# calc_tile_num<-function(lon_deg,lat_deg,zoom) {
#   ## Good notes: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Lon..2Flat._to_tile_numbers_2
#   n<-2^zoom
#   x<-n * ((lon_deg + 180)/360)
#   y_rad<-lat_deg * (pi/180)   #https://stackoverflow.com/questions/10140029/convert-latitude-longitude-in-degree-radians
#   y <- n * (1 - (log(tan(y_rad) + (1/cos(y_rad))) / pi)) / 2    #note that 1/cos is the secant function (secant function not available in base R)
#   xtile<-floor(x)
#   ytile<-floor(y)
#   ##-- keep the remainders -- useful for locating the point on the tile
#   xremainder<-x %% xtile
#   yremainder<-y %% ytile
#   return(c(xtile,ytile,xremainder,yremainder))
# }
#
# ####--- Find the corner points (in lat/lng) for the tile, to georeference the raster graphic.
# ## This is the xmin/ymin of the tile, plus the xmin/ymin of the next (+1) tile (to the right and down)
# find_tile_corners<-function(tile_x,tile_y,zoom) {
#   n<-2^zoom
#   #Tile origin
#   x_this_tile = tile_x / n * 360.0 - 180.0
#   y_this_tile_rad = atan(sinh(pi * (1 - 2 * tile_y / n)))
#   y_this_tile = y_this_tile_rad * 180.0 / pi
#   rm(y_this_tile_rad)
#   #Origin for over 1 tile, down 1 tile (+1 tile in each direction)
#   x_next_tile = (tile_x+1) / n * 360.0 - 180.0
#   y_next_tile_rad = atan(sinh(pi * (1 - 2 * (tile_y+1) / n)))
#   y_next_tile = y_next_tile_rad * 180.0 / pi
#   return(c(x_this_tile,x_next_tile,y_next_tile,y_this_tile))
# }
#
# ############ Form the API Call###########
# ## NOTE: Optional argument here to specify dates -- the default argument is Jan 1 2000 to Jan 1 2021; NM will return most recent image in that window.
# ## The argument tile should be a length-two vector with x and y tile numbers
# nm_api_call<-function(tile,zoom,apikey,
#                       time_start='2000-01-01',time_end='2021-01-01',mosaic='latest') {
#   #https://docs.nearmap.com/display/ND/Tile+API
#   call<-paste0('https://api.nearmap.com/tiles/v3/Vert/',
#                zoom,'/',tile[1],'/',tile[2],'.png?apikey=',apikey,
#                '&since=',time_start,'&until=',time_end,'&mosaic=',mosaic,'&tertiary=satellite')
#   return(call)
# }
#
# gettile<-function(tile_api_call,saveloc) {
#   curl::curl_download(tile_api_call,saveloc)  # can save it somewhere better than 'dest' in the working folder
#   a<-png::readPNG(saveloc)
#   return(a)
# }
#
# ############# Make the map figure in ggplot() ############
# # Use the png as a base_map with annotation_raster
# ## Helpful example: https://yutani.rbind.io/post/2018-06-09-plot-osm-tiles/
# ## Note for some reason, can't call ggplot() with no data -- have to put something in initial ggplot call.
#
# make_nearmap_figure<-function(baseimage,
#                               parcelsandthelike,
#                               tiles_corners,
#                               rasterimages,
#                               parcels,
#                               bldgpolygons,
#                               Final_Map_Points,
#                               roofs) {
#
#   ##-- find the plot limits based on the tiles
#   min_xmin<-min(tiles_corners[,'xmin'])
#   max_xmax<-max(tiles_corners[,'xmax'])
#   min_ymin<-min(tiles_corners[,'ymin'])
#   max_ymax<-max(tiles_corners[,'ymax'])
#
#   ##-- Initialize the map and choose pre-, post-imagery.
#   b<-ggplot(data=Final_Map_Points)
#   ##-- add each of the 4 tiles as a georeferenced raster image using corner points
#   if (baseimage==TRUE) {
#     for (i in 1:4) {
#       b<-b+annotation_raster(rasterimages[[i]],tiles_corners[i,1],tiles_corners[i,2],
#                              tiles_corners[i,3],tiles_corners[i,4])
#     }
#   }
#   ##-- add other elements to map
#   if (parcelsandthelike==TRUE) {
#     b <- b + geom_sf(data=parcels,color='gray',size=0.8,fill=NA) +
#       geom_sf(data=bldgpolygons,color='purple',fill=NA,alpha=0.5,size=0.3) +
#       geom_sf(data=roofs,color='blue',fill=NA,alpha=0.5,size=0.3) +
#       geom_sf(data=Final_Map_Points,aes(fill=as.character(destroyed),shape=as.character(destroyed)),color='white',stroke=1,size=1.25,alpha=1) +
#       scale_fill_viridis_d(option = "viridis",begin=0.5,breaks=c("TRUE","FALSE"))+
#       scale_shape_manual(values=c(21,22),breaks=c("TRUE","FALSE"))
#       #scale_color_manual(values=c('red','green'),breaks=c("TRUE","FALSE"))
#   }
#   ##-- Add last options to map
#   b <- b + theme_map() + theme(legend.position = "none") +
#     coord_sf(xlim = c(min_xmin,max_xmax), ylim = c(min_ymin,max_ymax), expand = FALSE)
#   return(b)
# }
#
# ############# Make a map figure with distances labeled ############
# make_nearmap_with_distances<-function(baseimage,
#                               parcelsandthelike,
#                               tiles_corners,
#                               rasterimages,
#                               parcels,
#                               bldgpolygons,
#                               Final_Map_Points,
#                               roofs,
#                               distmeasure='walltowall') {
#
#   ##-- find the plot limits based on the tiles
#   min_xmin<-min(tiles_corners[,'xmin'])
#   max_xmax<-max(tiles_corners[,'xmax'])
#   min_ymin<-min(tiles_corners[,'ymin'])
#   max_ymax<-max(tiles_corners[,'ymax'])
#
#   ##-- Initialize the map and choose pre-, post-imagery.
#   b<-ggplot(data=Final_Map_Points)
#   ##-- add each of the 4 tiles as a georeferenced raster image using corner points
#   if (baseimage==TRUE) {
#     for (i in 1:4) {
#       b<-b+annotation_raster(rasterimages[[i]],tiles_corners[i,1],tiles_corners[i,2],
#                              tiles_corners[i,3],tiles_corners[i,4])
#     }
#   }
#   ##-- add other elements to map
#   if (parcelsandthelike==TRUE) {
#     b <- b +
#       #geom_sf(data=parcels,color='gray',size=0.8,fill=NA) +
#       geom_sf(data=bldgpolygons,color='blue',fill=NA,alpha=0.5,size=0.3) +
#       geom_sf(data=Final_Map_Points,aes(shape=destroyed),size=1.5,alpha=0.5)
#   }
#   ##-- Add DISTANCE LABELS and other last options to map
#   if (distmeasure=='roofcentroid') {
#     b <- b + geom_sf_text(data=Final_Map_Points,aes(label=distfrom1),size=4.5,color='white') +
#       theme_map() + theme(legend.position = "none") +
#       coord_sf(xlim = c(min_xmin,max_xmax), ylim = c(min_ymin,max_ymax), expand = FALSE)
#   } else if (distmeasure=='walltowall') {
#     b <- b + geom_sf_text(data=roofs,aes(label=distfrom1),size=4.5,color='white') +
#       theme_map() + theme(legend.position = "none") +
#       coord_sf(xlim = c(min_xmin,max_xmax), ylim = c(min_ymin,max_ymax), expand = FALSE)
#   }
#   return(b)
# }
#
#
# #------------------------------------------------------------------#
# #--- Master function to run everything and save the figure --------#
# #------------------------------------------------------------------#
# ## Keep the nearmap option set to FALSE until ready to start downloading tiles - remember API download limits.
# # Prototype without the base maps to find a place that looks good for checking, then turn nearmap option to TRUE.
#
# qc_map<-function(importpcid,incid,nearmap=FALSE,zoom=18) {
#
#   ptm<-proc.time()
#   print("Setting up data and selecting building footprint polygons")
#
#   zerohome<-homes %>%
#     dplyr::filter(ImportParcelID==importpcid & incidentid==incid) %>%
#     dplyr::select(lon,lat,FIPS,ImportParcelID,incidentid,State)
#
#   ## Get the bounding box around the example home.
#
#   zoombox <- zerohome %>%
#     st_buffer(250) %>%   #This projection has units of meters, so this buffer should be measured in meters
#     st_bbox()
#
#   Final_Map_Points<-homes %>%
#     dplyr::filter(incidentid==incid) %>%  ## rare cases of homes inside zip codes affected by 2 WF incidents.
#     st_crop(zoombox)
#
#   ## -- calculate centroid distances
#   Final_Map_Points$distfrom1<-round(st_distance(Final_Map_Points[Final_Map_Points$ImportParcelID==importpcid,],
#                                                 Final_Map_Points,by_element=TRUE),0)
#
#
#   ##--- All building footprints that appear in the Microsoft imagery
#   bldgpolygons<-qread(file.path(WORKING,'buildings',paste0(incid,'bldgs.qs'))) %>%
#     st_transform(crs=st_crs(homes)) %>%
#     #st_make_valid() %>%
#     st_join(Final_Map_Points,left=FALSE)
#
#   ## Building Footprints identified as main residences and assigned locs in analysis
#   roofs<-roofshapes %>%
#     #st_make_valid() %>%
#     st_crop(zoombox)
#
#   ## -- Calculate wall to wall distances for these buildings (if I have a building footprint for the "zero" house)
#   if (nrow(roofs[roofs$ImportParcelID==importpcid,])>0) {
#     roofs$distfrom1<-round(st_distance(roofs[roofs$ImportParcelID==importpcid,],
#                                                 roofs,by_element=TRUE),0)
#     roofs$distfrom1[roofs$roof_fail==TRUE]<-NA
#   }
#
#   ### Parcel boundary lines
#   parcels<-lotlines %>%
#     st_crop(zoombox)
#
#   ### Fire perimeter
#   thisperim<-perims %>%
#     dplyr::filter(incidentid==incid) %>%
#     st_crop(zoombox)
#
#
#   ###### Put everything into WGS 84 for plotting
#   Final_Map_Points<-st_transform(Final_Map_Points,4326)
#   bldgpolygons<-st_transform(bldgpolygons,4326)
#   roofs<-st_transform(roofs,4326)
#   parcels<-st_transform(parcels,4326)
#   thisperim<-st_transform(thisperim,4326)
#
#   ###### Make a no-image schematic map to test
#   q<-ggplot() +
#     geom_sf(data=parcels,color='gray',fill=NA,alpha=0.5)+
#     geom_sf(data=bldgpolygons,color='yellow',fill=NA,alpha=0.5) +
#     geom_sf(data=roofs,color='blue',fill=NA,alpha=0.5) +
#     geom_sf(data=Final_Map_Points,aes(color=destroyed),size=1,alpha=0.5) +
#     geom_sf(data=thisperim,color='red',fill=NA) +
#     #geom_sf_text(data=Final_Map_Points,aes(label=ImportParcelID))+
#     theme_map()
#   save_plot(file.path(RES,'niceNMfigs',paste0(importpcid,'_',incid,'_lines.png')), plot = q)
#
#   if (nearmap==TRUE) {
#     print("Getting tile information and pinging NearMap API")
#
#     #### Get info for the tile that contains the point of interest
#     tileinfo<-calc_tile_num(zerohome[[1]],zerohome[[2]],zoom)
#     tile<-tileinfo[1:2]
#     ## Get info for the 3 adjacent tiles, on whichever side is closest to point
#     if (tileinfo[3]>0.5){
#       xmove <- 1
#     } else {
#       xmove<- -1
#     }
#     if (tileinfo[4]>0.5){
#       ymove <- 1
#     } else {
#       ymove<- -1
#     }
#     tiles<-list(tile,
#                 c(tile[1]+xmove,tile[2]),
#                 c(tile[1]+xmove,tile[2]+ymove),
#                 c(tile[1],tile[2]+ymove))
#
#     tiles_corners<-list(find_tile_corners(tiles[[1]][1],tiles[[1]][2],zoom),
#                         find_tile_corners(tiles[[2]][1],tiles[[2]][2],zoom),
#                         find_tile_corners(tiles[[3]][1],tiles[[3]][2],zoom),
#                         find_tile_corners(tiles[[4]][1],tiles[[4]][2],zoom))
#     tiles_corners<-data.frame(xmin=c(tiles_corners[[1]][1],tiles_corners[[2]][1],tiles_corners[[3]][1],tiles_corners[[4]][1]),
#                               xmax=c(tiles_corners[[1]][2],tiles_corners[[2]][2],tiles_corners[[3]][2],tiles_corners[[4]][2]),
#                               ymin=c(tiles_corners[[1]][3],tiles_corners[[2]][3],tiles_corners[[3]][3],tiles_corners[[4]][3]),
#                               ymax=c(tiles_corners[[1]][4],tiles_corners[[2]][4],tiles_corners[[3]][4],tiles_corners[[4]][4]))
#
#     #################################################################################
#     ############## PING NEARMAP API #################################################
#     #################################################################################
#
#     ### Load fire info
#     inc<-zerohome$incidentid
#     fire_ids<-readRDS(file.path(WORKING,"fire_ids.RDS"))
#
#     ####-- Prefire images
#     timeperiod_start<-'2000-01-01'
#     timeperiod_end<-as.character(fire_ids$StartDate[fire_ids$incidentid==inc])
#
#     api_calls<-lapply(tiles,nm_api_call,
#                       zoom=zoom,apikey=apikey,time_start=timeperiod_start,time_end=timeperiod_end,mosaic='latest') #mosaic = latest prioritizes most recent image within specified time window
#     print(api_calls)
#     preimages<-lapply(api_calls,gettile,saveloc=file.path(WORKING,"intermediatefiles","nmimage"))
#
#     ####-- Postfire images (unrestricted follow up period)
#     timeperiod_start<-as.character(fire_ids$StartDate[fire_ids$incidentid==inc])
#     timeperiod_end<-'2025-01-01'
#
#     api_calls<-lapply(tiles,nm_api_call,
#                       zoom=zoom,apikey=apikey,time_start=timeperiod_start,time_end=timeperiod_end,mosaic='earliest') #mosaic=earliest will preference earliest image within the specified time window (soonest after fire, in this case)
#     postimages<-lapply(api_calls,gettile,saveloc=file.path(WORKING,"intermediatefiles","nmimage"))
#
#     ####-- Restricted window Postfire images -- require within 365 days of fire; for assessing damage report accuracy
#     timeperiod_start<-as.character(fire_ids$StartDate[fire_ids$incidentid==inc])
#     timeperiod_end<-as.character((1*365)+fire_ids$StartDate[fire_ids$incidentid==inc])
#
#     api_calls<-lapply(tiles,nm_api_call,
#                       zoom=zoom,apikey=apikey,time_start=timeperiod_start,time_end=timeperiod_end,mosaic='earliest') #mosaic=earliest will preference earliest image within the specified time window (soonest after fire, in this case)
#     postimages_restricted<-lapply(api_calls,gettile,saveloc=file.path(WORKING,"intermediatefiles","nmimage"))
#
#
#     #################################################################################
#     ############## MAKE THE FIGURE #################################################
#     #################################################################################
#
#     print("Making and saving figures")
#
#     namestem<-importpcid
#
#     #-- "before" imagery
#     pre<-make_nearmap_figure(baseimage=TRUE,
#                              parcelsandthelike=TRUE,
#                              tiles_corners,
#                              preimages,
#                              parcels,
#                              bldgpolygons,
#                              Final_Map_Points,
#                              roofs)
#     save_plot(file.path(RES,"niceNMfigs",paste0(namestem,"_pre.png")), plot = pre, base_asp =1)
#
#     #-- "after" imagery
#     post<-make_nearmap_figure(baseimage=TRUE,
#                               parcelsandthelike=TRUE,
#                               tiles_corners,
#                               postimages,
#                               parcels,
#                               bldgpolygons,
#                               Final_Map_Points,
#                               roofs)
#     save_plot(file.path(RES,"niceNMfigs",paste0(namestem,"_post.png")), plot = post, base_asp =1)
#
#     #-- "after" imagery with restricted time window
#     post<-make_nearmap_figure(baseimage=TRUE,
#                               parcelsandthelike=TRUE,
#                               tiles_corners,
#                               postimages_restricted,
#                               parcels,
#                               bldgpolygons,
#                               Final_Map_Points,
#                               roofs)
#     save_plot(file.path(RES,"niceNMfigs",paste0(namestem,"_post_restricted.png")), plot = post, base_asp =1)
#
#     #-- no imagery
#     noimage<-make_nearmap_figure(baseimage=FALSE,
#                                  parcelsandthelike=TRUE,
#                                  tiles_corners,
#                                  preimages,
#                                  parcels,
#                                  bldgpolygons,
#                                  Final_Map_Points,
#                                  roofs)
#     save_plot(file.path(RES,"niceNMfigs",paste0(namestem,"_noimage.png")), plot = noimage, base_asp =1)
#     #-- no parcel lines etc
#     nolines<-make_nearmap_figure(baseimage=TRUE,
#                                  parcelsandthelike=FALSE,
#                                  tiles_corners,
#                                  postimages,
#                                  parcels,
#                                  bldgpolygons,
#                                  Final_Map_Points,
#                                  roofs)
#     save_plot(file.path(RES,"niceNMfigs",paste0(namestem,"_noparcels.png")), plot = nolines, base_asp =1)
#
#     #-- Map with distances to neighbors: Roof centroids
#     withdists<-make_nearmap_with_distances(baseimage=TRUE,
#                                            parcelsandthelike=TRUE,
#                                            tiles_corners,
#                                            preimages,
#                                            parcels,
#                                            bldgpolygons,
#                                            Final_Map_Points,
#                                            roofs,
#                                            distmeasure='roofcentroid')
#     save_plot(file.path(RES,"niceNMfigs",paste0(namestem,"_dist_centroids.png")), plot = withdists, base_asp =1)
#
#
#     #-- Map with distances to neighbors: Wall to wall (if I have a footprint for the "zero" house)
#     if (nrow(roofs[roofs$ImportParcelID==importpcid,])>0) {
#       withdists<-make_nearmap_with_distances(baseimage=TRUE,
#                                              parcelsandthelike=TRUE,
#                                              tiles_corners,
#                                              preimages,
#                                              parcels,
#                                              bldgpolygons,
#                                              Final_Map_Points,
#                                              roofs,
#                                              distmeasure='walltowall')
#       save_plot(file.path(RES,"niceNMfigs",paste0(namestem,"_dist_walls.png")), plot = withdists, base_asp =1)
#     }
#
#   }
#   return(paste0('Run time ',(proc.time()-ptm)[[3]] /60, ' minutes'))
# }
#
#
#
#
#
# #---------------------------------------------------------#
# #------- Get NearMap Coverage Polygons -------------------#
# #---------------------------------------------------------#
#
# ## Useful for selecting homes that will have imagery
#
# z<-homes %>%
#   st_transform(crs=4326) %>%
#   st_bbox() %>%
#   st_as_sfc() %>%
#   st_coordinates() %>%
#   data.frame()
#
#
# zz<-character()
# for (i in 1:nrow(z)){
#   zz<-c(zz,z$X[i],z$Y[i])
# }
# zz<-paste(zz,collapse=",")
#
# covgapicall<-paste0(
#   'https://api.nearmap.com/coverage/v2/surveyresources/boundaries.geojson?&apikey=',
#   apikey,
#   '&polygon=',
#   zz,
#   '&limit=10000&resources=tiles:Vert')
#
# curl::curl_download(covgapicall,destfile=file.path(WORKING,'nmcovg.json'))
#
# covg<-st_read(file.path(WORKING,'nmcovg.json'))
#
# #ggplot(data=st_transform(homes,crs=st_crs(covg)))+
# #         geom_sf(aes(color=incidentid),alpha=0.25)+
# #         geom_sf(data=covg,color='blue',fill=NA,alpha=0.33)+
# #         theme_map()+theme(legend.position='none')
#
#
# #------------------------------------------------------------------#
# #------- Check data usage and rate limit status -------------------#
# #------------------------------------------------------------------#
#
# #-- Show API headers, including rate limit, current usage, and reset date
# f<-curlGetHeaders(paste0('https://api.nearmap.com/tiles/v3/Vert/19/119799/215845.jpg?apikey=',apikey))
# f
#
# #-- Convert Reset time from UNIX epoch to human date
# reset_epoch<-as.numeric(gsub('[^0-9]','',f[14]))
# as.POSIXct(reset_epoch, origin="1970-01-01")
#
#
# #------------------------------------------------------------------#
# #------ Examples used in paper ------------------------------------#
# #------------------------------------------------------------------#
#
# carr1<-c('153102589','CASHU 007808')
# carr2<-c('153102455','CASHU 007808')
# lilac1<-c('17684182','CAMVU 024612')
# tubbs1<-c('20577669','CALNU 010104')
# thomas1<-c('21205528','CAVNC 103156')
#
# toplot<-list(carr1,carr2,lilac1,tubbs1,thomas1)
#
#
# ### --- Fig XX: Example distances
# for (i in 1:length(toplot)) {
#   qc_map(importpcid = toplot[[i]][1],
#          incid=toplot[[i]][2],
#          nearmap=TRUE,zoom=18)
# }
#
# ## Copy the files actually used in paper to main results folder
# forpaper<-c('153102589_pre.png','153102589_post.png',
#             '20577669_dist_walls.png','21205528_dist_walls')
# a<-list.files(file.path(RES,'niceNMfigs'),
#                     pattern=paste(forpaper,collapse='|'))
# file.rename(from=file.path(RES,'niceNMfigs',a),
#             to=file.path(RES,a))
#
# #--------------------------------------------------------------------------------#
# #### For a given incident, map a handful of homes at random
# #--------------------------------------------------------------------------------#
#
# m<-'CACZU 005205'
#
# foo<-homes %>%
#   filter(incidentid==m) %>%
#   st_transform(st_crs(covg)) %>%
#   st_join(covg,left=FALSE) %>%
#   select(ImportParcelID) %>%
#   slice_sample(n=8)
#
# for (k in 1:nrow(foo)) {
#   print(k)
#   qc_map(importpcid = foo$ImportParcelID[k],
#          incid=m,
#          nearmap=TRUE,zoom=18)
# }
#
#
#
#
# #--------------------------------------------------------------------------------#
# #### Quality Checking Locations and Damage
# #--------------------------------------------------------------------------------#
#
#
# # TODO: NOT SURE IF THIS MATCHES THE CURRENT VARIABLE NAMES
# density<-read_fst(file.path(WORKING, "neighbor-vars.fst")) %>%
#   select(ImportParcelID,UNIT_AND_NUMBER,BuildingOrImprovementNumber,starts_with('totwithin')) %>%
#   rename(incidentid=UNIT_AND_NUMBER)
#
# loc_sources<-readRDS(file.path(WORKING,"zillowdata.RDS")) %>%
#   select(ImportParcelID,loc_source) %>%
#   distinct(ImportParcelID,.keep_all=T)
#
# modestatistic<-function(x) which.max(tabulate(x))
#
# set.seed(3285972)
#
# foo<-homes %>%
#   filter(State=="California") %>%
#   select(ImportParcelID,incidentid,BuildingOrImprovementNumber) %>%
#   st_transform(st_crs(covg)) %>%
#   st_join(covg,left=FALSE) %>%
#   st_drop_geometry() %>%
#   distinct(ImportParcelID,incidentid,BuildingOrImprovementNumber) %>%  ## Merging to coverage tiles creates duplicates (multiple coverage polygons for 1 house)
#   merge(density,by=c('ImportParcelID','BuildingOrImprovementNumber','incidentid')) %>%
#   filter(!is.na(totwithin200m)) %>% ## Missing on this variable means not usable in neighbor analysis (e.g., condo)
#   mutate(bin=cut(totwithin200m,breaks=c(0,10,15),include.lowest=TRUE)) %>%
#   merge(loc_sources,by='ImportParcelID') %>%
#   mutate(loc_source=factor(loc_source)) %>%
#   ## Identify the main location source for each incident
#   group_by(incidentid) %>%
#   mutate(dominantlocmethod=modestatistic(loc_source),
#          dominantlocmethod=levels(loc_source)[dominantlocmethod]) %>%
#   filter(loc_source==dominantlocmethod) %>%
#   group_by(incidentid,bin,dominantlocmethod) %>%
#   slice_sample(n=1)
#
# dir.create(file.path(RES, "niceNMfigs"), showWarnings = F)
# resave<-function(f) {
#   fname<-list.files(file.path(RES,'niceNMfigs'),pattern=paste0(foo$ImportParcelID[k],f))
#   file.rename(from=file.path(RES,'niceNMfigs',fname),to=file.path(RES,'qc','roundN',fname))
# }
#
# dir.create(file.path(RES, "qc", "roundN"), showWarnings = F)
# datasheet<-foo %>%
#   select(ImportParcelID,incidentid,totwithin200m,bin,dominantlocmethod) %>%
#   write.csv(file.path(RES,'qc','roundN','datasheet.csv'))
#
# for (k in c(1:80,82:nrow(foo))) {  #I don't know why image 81 and 82 doesn't work? Gives 404 error.
# #for (k in 81:82) {
#   print(k)
#   try(
#     qc_map(importpcid = foo$ImportParcelID[k],
#          incid=foo$incidentid[k],
#          nearmap=TRUE,zoom=18)
#     ,TRUE)
#   ## Resave needed images elsewhere and cleanup the NM results folder
#   resave('_pre')
#   resave('_post')
#   resave('_post_restricted')
#   unlink(list.files(file.path(RES,'niceNMfigs'),pattern=as.character(foo$ImportParcelID[k]),full.names = TRUE))
# }

#### -- Read in the hand-evaluated accuracy scores and make a table ----

tempfunc <- function(x) ifelse(x == "na", NA, x)

data <- read.csv(file.path(INPUT_PRIVATE, "qc", "datasheet_handevaluated.csv")) %>%
  mutate(across(matches("Total|BadLoc|BadDamage"), ~ as.numeric(tempfunc(.x))))

fire_ids <- readRDS(file.path(WORKING, "fire_ids.RDS"))

d <- merge(data, fire_ids, by = "incidentid", all.x = TRUE) %>%
  mutate(
    year = year(StartDate),
    LocationErrorRate = BadLoc / Total,
    DamageErrorRate = BadDamage / Total
  )

loc <- d %>%
  select(
    incidentid, ImportParcelID, bin, dominantlocmethod, Total,
    LocationErrorRate, HasPre, year
  ) %>%
  filter(HasPre == "yes") %>%
  group_by(bin, dominantlocmethod) %>%
  summarize(
    LocationError = round(weighted.mean(LocationErrorRate, w = Total, na.rm = TRUE), 3),
    images = n(),
    count = sum(Total, na.rm = TRUE)
  ) %>%
  ungroup() %>%
  arrange(dominantlocmethod, bin) %>%
  mutate(group = case_when(
    bin == "(10,15]" & dominantlocmethod == "BldgFootprints" ~ "Footprint Geocodes, 10+ within 200 m",
    bin == "[0,10]" & dominantlocmethod == "BldgFootprints" ~ "Footprint Geocodes, 0--9 within 200 m",
    bin == "(10,15]" & dominantlocmethod == "ESRI" ~ "Address Geocodes, 10+ within 200 m",
    bin == "[0,10]" & dominantlocmethod == "ESRI" ~ "Address Geocodes, 0--9 within 200 m"
  )) %>%
  select(group, LocationError, count, images) %>%
  transmute(` ` = group, Images = images, Homes = count, `Error Rate` = LocationError * 100)

dam <- d %>%
  select(incidentid, ImportParcelID, Total, DamageErrorRate, HasPostRest, year) %>%
  ungroup() %>%
  summarize(
    DamageError = round(weighted.mean(DamageErrorRate, w = Total, na.rm = TRUE), 3),
    images = n(),
    count = sum(Total, na.rm = TRUE)
  ) %>%
  mutate(group = "All") %>%
  select(group, DamageError, count, images) %>%
  transmute(` ` = group, Images = images, Homes = count, `Error Rate` = DamageError * 100)

error_rates_df <- bind_rows(loc, dam) %>%
  transmute(` `, `Error (\\%)` = `Error Rate`, Homes, Images)

kbl(error_rates_df,
    format = "latex",
    escape = F,
    align = "lrrr",
    booktabs = T,
    format.args = list(big.mark = ","),
    linesep = "",
    digits = 3
) %>%
  group_rows("Panel A. Location Error Rates", 1, 4, bold = F, italic = T) %>%
  group_rows("Panel B. Damage Error Rate", 5, 5, bold = F, italic = T) %>%
  save_kable(file = "results/error-rates.tex")



##### California Examples for Slides, Etc #####################


# qc_map("11592188",
#   incid = "CALAC 326166",
#   nearmap = TRUE, zoom = 18
# )
# 
# 
# morehomes <- data.frame(
#   homes = c(
#     "21189269",
#     "20704385",
#     "11171486",
#     "153112444",
#     "153063653"
#   ),
#   fires = c(
#     "CAVNC 103156",
#     "CALNU 010104",
#     "CALAC 33891",
#     "CASHU 007808",
#     "CASHU 007808"
#   )
# )
# 
# for (k in 1:nrow(morehomes)) { # nrow(morehomes)
#   qc_map(
#     importpcid = morehomes$homes[k],
#     incid = morehomes$fires[k],
#     nearmap = TRUE, zoom = 18
#   )
# }
# 
# ### ALREADY RUN EXAMPLES:
# 
# # Shasta County example
# qc_map(
#   importpcid = "153102118",
#   incid = "CASHU 007808",
#   nearmap = TRUE, zoom = 18
# )
# 
# # Better Ventura County example: Thomas Fire
# qc_map("156545852",
#   incid = "CAVNC 103156",
#   nearmap = TRUE, zoom = 18
# )
# 
# # Sonoma example (Santa Rosa)
# qc_map("20577669",
#   incid = "CALNU 010104",
#   nearmap = TRUE, zoom = 18
# )
# 
# # Ventura County example
# qc_map("156545852",
#   incid = "CAVNC 103156",
#   nearmap = TRUE, zoom = 18
# )
# 
# # Woolsey Fire example
# qc_map("151081507",
#   incid = "CALAC 33891",
#   nearmap = TRUE, zoom = 18
# )
# 
# 
# # San Diego example
# importpcid <- "18659806"
# incid <- "Witch2007"
# nearmap <- TRUE
# zoom <- 18
# 
# qc_map("11592188",
#   incid = "CALAC 326166",
#   nearmap = TRUE, zoom = 18
# )
